################################################################################
## REPLICATION R CODE FOR                                                     ##
##                                                                            ##
## ROD ALENCE AND XICHAVO ALECIA NDLOVU,                                      ##
## "POLITICAL ACCOUNTABILITY AND DEVELOPMENT IN AFRICA'S RESOURCE ECONOMIES", ##
## *THE EXTRACTIVE INDUSTRIES AND SOCIETY* 22 (2025): 101634.                 ##
################################################################################

#############
## SESSION ##
#############
suppressPackageStartupMessages({
library(wordcloud)
library(stargazer)
library(marginaleffects)
library(systemfit)
library(lmtest)
library(cv)
library(dplyr)
library(ggplot2)
library(WhatIf)
})
load("exis_2025.RData")
sessionInfo()

####################
## TRANSFORM DATA ##
####################
## Logarithmic transformations
exis$lgdpnrPC9095 <- log(exis$gdpnrPC9095)
exis$lrentPC9095 <- log(exis$rentPC9095 + 1)
exis$lgdpPC9095 <- log(exis$gdpPC9095)

## Standardize variables
Stdize <- function(x) res <- (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)

exis$Z.lgdpnrPC9095 <- Stdize(exis$lgdpnrPC9095)
exis$Z.lrentPC9095 <- Stdize(exis$lrentPC9095)
exis$Z.incl9095 <- Stdize(exis$incl9095)
exis$Z.incl1519 <- Stdize(exis$incl1519)
exis$Z.sust95 <- Stdize(exis$sust95)
exis$Z.sust1518 <- Stdize(exis$sust1518)
exis$Z.dem_avg <- Stdize(exis$dem_avg)
exis$Z.compete_avg <- Stdize(exis$compete_avg)
exis$Z.parlinks_avg <- Stdize(exis$parlinks_avg)
exis$Z.war_avg <- Stdize(exis$war_avg)

###################
## PLOT FIGURE 2 ##
###################
## FIGURE 2: DEVELOPMENT TRAJECTORIES
exis$incl_change <- exis$incl1519 - exis$incl9095
exis$sust_change <- exis$sust1518 - exis$sust95
## Bivariate correlation:
cor(exis[, c("incl_change", "sust_change")], use = "complete.obs") |> round(2)

## Plot figure
par(mar = c(3.1, 3.1, 0.6, 0.6),
    cex.axis = 0.7,
    cex.lab = 0.9,
    mgp = c(2, 0.5, 0))
with(exis[!is.na(exis$sust_change), ], {
  plot(incl_change, sust_change, type = "n",
       xlab = "Social inclusiveness (change, 1990-95 to 2015-19)",
       ylab = "Economic sustainability (change, 1995 to 2015-18)",
       xlim = c(0, 0.35),
       ylim = c(-0.3, 0.95),
       axes = FALSE)
  axis(side = 1, at = seq(0, 0.3, by = 0.1))
  axis(side = 2, at = seq(-0.3, 0.9, by = 0.3))
  abline(lm(sust_change ~ incl_change), lty = "dotted")
  mtext("Pearson's r = 0.24", line = -1, at = 0.05, cex = 0.9)
  textplot(incl_change, sust_change, words = country, cex = 0.75,
           show.lines = FALSE, new = FALSE)
})

###################
## PLOT FIGURE 3 ##
###################
## Party-citizen links (rank)
exis$parlinks21_rank <- nrow(exis) - rank(exis$parlinks21, na.last = "keep") + 1
## Country labels that reflect "resource intensity" in all-caps
exis$resIMF21_lab <- exis$country
exis[exis$resIMF21, "resIMF21_lab"] <- toupper(exis[exis$resIMF21, "resIMF21_lab"])
## Regime category labels
regime_labs <- c("",
                 "Non-democracy",
                 "Dominant-party\ndemocracy",
                 "Competitive\ndemocracy",
                 "")
par(mar = c(4.1, 3.1, 1.1, 0.1), cex.axis = 0.8)
plot(as.numeric(exis$polregime21), exis$parlinks21_rank,
     xlim = c(0.4, 3.5),
     ylim = c(47, 1),
     type="n",
     axes = FALSE,
     ylab = "",
     xlab = "Political regime, 2021")
segments(x0=c(0.4, 1.2), x1 = c(0.8, 3.5), y0=c(24, 24),
         lwd = 1, lty = "16")
text(exis$polregime21, exis$parlinks21_rank, labels = exis$resIMF21_lab,
     cex = 0.7)
axis(side = 2, at = c(1, 12, 24, 36, 47),
     lwd = 1, col = "grey60", lwd.ticks = 1,
     mgp = c(0.5, 0.5, 0))
axis(side = 1, at = c(0.4, 1:3, 3.5),
     labels = regime_labs,
     mgp = c(3, 0.5, 0), padj = 0.5,
     lwd = 1, col = "grey60", lwd.ticks = 1)
mtext("Programmatic party linkages, 2021 (rank)",
      side = 2, line = 2)
mtext('"Resource-intensive" economies (IMF) in capitals',
      side = 3, at = 0.4, line = 0, adj = 0, cex = 0.8, font = 3)

#################################
## RUN REGRESSIONS FOR TABLE 1 ##
#################################
## Inclusiveness as outcome (columns 1-3)
exis_incl <- exis
exis_incl$lagged_outcome <- exis_incl$Z.incl9095  # for stargazer

incl_1 <- lm(Z.incl1519 ~ Z.lgdpnrPC9095 + Z.lrentPC9095,
             data = exis_incl)
summary(incl_1)

incl_2 <- lm(Z.incl1519 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095,
             data = exis_incl)
summary(incl_2)

incl_3 <- lm(Z.incl1519 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg,
             data = exis_incl)
summary(incl_3)

## Sustainability as outcome (columns 4-6)
exis_sust <- exis
exis_sust$lagged_outcome <- exis_incl$Z.sust95  # for stargazer

sust_1 <- lm(Z.sust1518 ~ Z.lgdpnrPC9095 + Z.lrentPC9095,
             data = exis_sust)
summary(sust_1)

sust_2 <- lm(Z.sust1518 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095,
             data = exis_sust)
summary(sust_2)

sust_3 <- lm(Z.sust1518 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg,
             data = exis_sust)
summary(sust_3)

## Regression table (LaTeX)
cov_labs <- c("Lagged outcome", "Non-rents/capita (log)",
              "Rents/capita (log)",
              "Electoral democracy", "Competitive", "Party linkages")
dv_labs <- c("Inclusiveness", "Sustainability")
table_note <- paste("OLS estimates and standard errors;",
                    "all variables standardized (z-scores).")

stargazer(incl_1, incl_2, incl_3, sust_1, sust_2, sust_3,
          type = "latex",
          title = "Predictors of development outcomes, 1990--2019",
          align = TRUE,
          keep.stat = c("n", "adj.rsq"),
          df = FALSE,
          covariate.labels = cov_labs,
          dep.var.labels = dv_labs,
          dep.var.caption = "Development outcomes",
          digits = 2,
          star.cutoffs = NA,
          omit.table.layout = "n")

##################################################
## RUN (AND SUMMARIZE) SIMULATIONS FOR FIGURE 4 ##
##################################################
R_n <- 100       # number of replicates per profile
## Random seeds for each dimension (to avoid inducing association)
incl_seed <- 53  # random seed -- inclusiveness
sust_seed <- 22  # random seed -- sustainability
## Function to extract draws
getDraws <- function(x) x[["draw"]]

## Setting hypothetical values to use for profiles
hyp_names <- c("low", "high")
hyp_tail <- 0.18

dem_hyp <- quantile(exis$Z.dem_avg, c(hyp_tail, 1 - hyp_tail), na.rm = TRUE)
names(dem_hyp) <- hyp_names

comp_hyp <- quantile(exis$Z.compete_avg, c(hyp_tail, 1 - hyp_tail), na.rm = TRUE)
names(comp_hyp) <- hyp_names

party_hyp <- quantile(exis$Z.parlinks_avg, c(hyp_tail, 1 - hyp_tail), na.rm = TRUE)
names(party_hyp) <- hyp_names

rent_hyp <- c(-0.31, 0.69)
names(rent_hyp) <- hyp_names

## Convert rent values back to dollars
invRents <- function(x) {
  res <- exp(mean(exis$lrentPC9095, na.rm = TRUE) +
               (x * sd(exis$lrentPC9095, na.rm = TRUE))) - 1
  return(res)
}
invRents(rent_hyp)

## Values for each political regime configuration (see Table 2)

## Competitive democracy with programmatic parties (upper-right)
dem_nd <- data.frame(
  lagged_outcome = 0,
  Z.lgdpnrPC9095 = 0,
  Z.dem_avg = dem_hyp["high"],
  Z.parlinks_avg = party_hyp["high"],
  Z.compete_avg = comp_hyp["high"]
)
## Dominant-party democracy (lower-right)
demDom_nd <- data.frame(
  lagged_outcome = 0,
  Z.lgdpnrPC9095 = 0,
  Z.dem_avg = dem_hyp["high"],
  Z.parlinks_avg = party_hyp["high"],
  Z.compete_avg = 0
)
## Competitive democracy (upper-left)
demComp_nd <- data.frame(
  lagged_outcome = 0,
  Z.lgdpnrPC9095 = 0,
  Z.dem_avg = dem_hyp["high"],
  Z.parlinks_avg = 0, 
  Z.compete_avg = comp_hyp["high"]
)
## Non-democracy (lower-left)
nondem_nd <- data.frame(
  lagged_outcome = 0,
  Z.lgdpnrPC9095 = 0,
  Z.dem_avg = dem_hyp["low"],
  Z.parlinks_avg = party_hyp["low"],
  Z.compete_avg = comp_hyp["low"]
)

## Modify each to get "resource-poor" (P) and "resource-rich" (R) variants
dem_P_nd <- dem_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["low"])
dem_R_nd <- dem_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["high"])

demDom_P_nd <- demDom_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["low"])
demDom_R_nd <- demDom_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["high"])

demComp_P_nd <- demComp_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["low"])
demComp_R_nd <- demComp_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["high"])

nondem_P_nd <- nondem_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["low"])
nondem_R_nd <- nondem_nd |>
  mutate(Z.lrentPC9095 = rent_hyp["high"])

## Run simulations (with marginaleffects::predictions)

## RESOURCE-POOR PROFILES:
## Competitive/programmatic
set.seed(incl_seed)
incl_dem_P_sims <- predictions(incl_3,
                             newdata = dem_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_dem_P_sims)

set.seed(sust_seed)
sust_dem_P_sims <- predictions(sust_3,
                              newdata = dem_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_dem_P_sims)

## Nondemocracy
set.seed(incl_seed)
incl_nondem_P_sims <- predictions(incl_3,
                             newdata = nondem_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_nondem_P_sims)

set.seed(sust_seed)
sust_nondem_P_sims <- predictions(sust_3,
                              newdata = nondem_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
set.seed(incl_seed)
summary(sust_nondem_P_sims)

## Competitive/clientelistic
incl_demComp_P_sims <- predictions(incl_3,
                             newdata = demComp_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_demComp_P_sims)

set.seed(sust_seed)
sust_demComp_P_sims <- predictions(sust_3,
                              newdata = demComp_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_demComp_P_sims)

## Dominant-programmatic
set.seed(incl_seed)
incl_demDom_P_sims <- predictions(incl_3,
                             newdata = demDom_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_demDom_P_sims)

set.seed(sust_seed)
sust_demDom_P_sims <- predictions(sust_3,
                              newdata = demDom_P_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_demDom_P_sims)

## RESOURCE-RICH PROFILES:
## Competitive/programmatic
set.seed(incl_seed)
incl_dem_R_sims <- predictions(incl_3,
                             newdata = dem_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_dem_R_sims)

set.seed(sust_seed)
sust_dem_R_sims <- predictions(sust_3,
                              newdata = dem_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_dem_R_sims)

## Nondemocracy
set.seed(incl_seed)
incl_nondem_R_sims <- predictions(incl_3,
                             newdata = nondem_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_nondem_R_sims)

set.seed(sust_seed)
sust_nondem_R_sims <- predictions(sust_3,
                              newdata = nondem_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_nondem_R_sims)

set.seed(incl_seed)

## Competitive/clientelistic
incl_demComp_R_sims <- predictions(incl_3,
                             newdata = demComp_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_demComp_R_sims)

set.seed(sust_seed)
sust_demComp_R_sims <- predictions(sust_3,
                              newdata = demComp_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_demComp_R_sims)

## Dominant-programmatic
set.seed(incl_seed)
incl_demDom_R_sims <- predictions(incl_3,
                             newdata = demDom_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(incl_demDom_R_sims)

set.seed(sust_seed)
sust_demDom_R_sims <- predictions(sust_3,
                              newdata = demDom_R_nd) |>
  inferences(method = "simulation", R = R_n) |>
  posterior_draws() |>
  getDraws()
summary(sust_demDom_R_sims)

###################
## PLOT FIGURE 4 ##
###################
bo_pal <- c("#4A5AB5", "B5A54A")  # color palette
## Faceted data frame with simulations
facet_polrent <- data.frame(
  type = factor(rep(1:4, each = R_n * 2),
                labels = c("Competitive Democracy/\nModerate Party Linkages",
                           "Competitive Democracy/\nProgrammatic Parties",
                           "Nondemocracy",
                           "Moderate Competitiveness/\nProgrammatic Parties")),
  rent = factor(rep(rep(1:2, each = R_n), 4),
                labels = c("Resource poor", "Resource rich")),
  incl = c(incl_demComp_P_sims, incl_demComp_R_sims,
           incl_dem_P_sims, incl_dem_R_sims,
           incl_nondem_P_sims, incl_nondem_R_sims,
           incl_demDom_P_sims, incl_demDom_R_sims),
  sust = c(sust_demComp_P_sims, sust_demComp_R_sims,
           sust_dem_P_sims, sust_dem_R_sims,
           sust_nondem_P_sims, sust_nondem_R_sims,
           sust_demDom_P_sims, sust_demDom_R_sims)
)

ggplot(data = facet_polrent, aes(x = incl, y = sust)) +
  geom_point(data = select(facet_polrent, -type), color= "grey80", size = 0.5) +
  geom_point(size=0.7, fill = "white",
             aes(shape = rent, color = rent)) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_vline(xintercept=0, linetype = "dashed") +
  scale_shape_manual(values=c(21, 3))+
  scale_color_manual(values =  c("Resource poor" = "#0072B2",
                                 "Resource rich" = "#D55E00")) +
  facet_wrap(~ type, nrow = 2) +
  theme_bw() +
  theme(legend.title = element_blank(),
        strip.background = element_rect(fill = "white"),
        strip.text = element_text(size = 9),
        legend.text = element_text(size = 9),
        legend.position = "top",
        legend.direction = "vertical",
        axis.text = element_text(size = 8)) +
  xlab("Social inclusiveness") +
  ylab("Economic sustainability")

#####################################################
## CHECK PLAUSIBILITY OF "COUNTERFACTUAL" PROFILES ##
#####################################################
## Non-democracy
(nondem_P_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                       data = exis, cfact = nondem_P_nd,
                       distance = "euclidian"))
(nondem_R_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                       data = exis, cfact = nondem_R_nd,
                       distance = "euclidian"))
## Competitive democracy
(demComp_P_cf <- whatif(~ Z.lrentPC9095 +
                          Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                        data = exis, cfact = demComp_P_nd,
                        distance = "euclidian"))
(demComp_R_cf <- whatif(~ Z.lrentPC9095 +
                          Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                        data = exis, cfact = demComp_R_nd,
                        distance = "euclidian"))
## Dominant-party democracy
(demDom_P_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                       data = exis, cfact = demDom_P_nd,
                       distance = "euclidian"))
(demDom_R_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                       data = exis, cfact = demDom_R_nd,
                       distance = "euclidian",
                       return.distance = TRUE))

## Full democracy
(dem_P_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                    data = exis, cfact = dem_P_nd,
                    distance = "euclidian"))
(dem_R_cf <- whatif(~ Z.lrentPC9095 + Z.dem_avg + Z.compete_avg +
                        Z.parlinks_avg,
                    data = exis, cfact = dem_R_nd,
                    distance = "euclidian",
                    return.distance = TRUE))

###############################################################
## REPORT ALTERNATIVE SPECIFICATION MENTIONED IN THE ARTICLE ##
###############################################################

#####################
## Robustness: War ##
#####################
## War-years as a predictor (4a) and "no-war" subset (4b)
## Inclusiveness
incl_4a <- lm(Z.incl1519 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg + Z.war_avg,
             data = exis_incl)
summary(incl_4a)  # include proportion war years as predictor

incl_4b <- lm(Z.incl1519 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg,
             data = subset(exis_incl, war_avg == 0))
summary(incl_4b)  # subset sample to countries with no civil war
## Sustainability
sust_4a <- lm(Z.sust1518 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg + Z.war_avg,
             data = exis_sust)
summary(sust_4a)  # include proportion war years as predictor

sust_4b <- lm(Z.sust1518 ~ lagged_outcome + Z.lgdpnrPC9095 + Z.lrentPC9095 +
            Z.dem_avg + Z.compete_avg + Z.parlinks_avg,
             data = subset(exis_incl, war_avg == 0))
summary(sust_4b)  # subset sample to countries with no civil war

#######################################
## Interactions and cross-validation ##
#######################################
## Outcome: Inclusiveness
## Additive
incl_add <- lm(Z.incl1519 ~  Z.incl9095 + Z.lgdpnrPC9095 + Z.lrentPC9095 +
                      (Z.dem_avg + Z.parlinks_avg + Z.compete_avg),
                data = exis)
summary(incl_add)
(incl_add_mse <- cv(incl_add, k = "loo")[["CV crit"]][1])

## Political interactions only
incl_intP <- lm(Z.incl1519 ~ Z.incl9095 + Z.lgdpnrPC9095 + Z.lrentPC9095 +
                       (Z.dem_avg * Z.parlinks_avg * Z.compete_avg),
                     data = exis)
summary(incl_intP)
(incl_intP_mse <- cv(incl_intP, k = "loo")[["CV crit"]][1])

## Rent interactions only
incl_intR <- lm(Z.incl1519 ~ Z.incl9095 + Z.lgdpnrPC9095 + Z.lrentPC9095 *
                       (Z.dem_avg + Z.compete_avg + Z.parlinks_avg),
                     data = exis)
summary(incl_intR)
(incl_intR_mse <- cv(incl_intR, k = "loo")[["CV crit"]][1])

## "All" interactions
incl_intA <- lm(Z.incl1519 ~ Z.incl9095 + Z.lgdpnrPC9095 + Z.lrentPC9095 *
                       (Z.dem_avg * Z.parlinks_avg * Z.compete_avg),
                     data = exis)
summary(incl_intA)
(incl_intA_mse <- cv(incl_intA, k = "loo")[["CV crit"]][1])

## Outcome: Sustainability
## Additive
sust_add <- lm(Z.sust1518 ~ Z.sust95 + Z.lgdpnrPC9095 + Z.lrentPC9095 +
                      (Z.dem_avg + Z.parlinks_avg + Z.compete_avg),
                data = exis)
summary(sust_add)
(sust_add_mse <- cv(sust_add, k = "loo")[["CV crit"]][1])

## Political interactions only
sust_intP <- lm(Z.sust1518 ~ Z.sust95 + Z.lgdpnrPC9095 + Z.lrentPC9095 +
                       (Z.dem_avg * Z.parlinks_avg * Z.compete_avg),
                     data = exis)
summary(sust_intP)
(sust_intP_mse <- cv(sust_intP, k = "loo")[["CV crit"]][1])

## Rent interactions only
sust_intR <- lm(Z.sust1518 ~ Z.sust95 + Z.lgdpnrPC9095 + Z.lrentPC9095 *
                       (Z.dem_avg + Z.compete_avg + Z.parlinks_avg),
                     data = exis)
summary(sust_intR)
(sust_intR_mse <- cv(sust_intR, k = "loo")[["CV crit"]][1])

## "All" interactions
sust_intA <- lm(Z.sust1518 ~ Z.sust95 + Z.lgdpnrPC9095 + Z.lrentPC9095 *
                       (Z.dem_avg * Z.parlinks_avg * Z.compete_avg),
                     data = exis)
summary(sust_intA)
(sust_intA_mse <- cv(sust_intA, k = "loo")[["CV crit"]][1])

####################################
## Seemingly unrelated regression ##
####################################
incl <- incl1519 ~ incl9095 + lgdpnrPC9095 + lrentPC9095 + dem_avg +
  parlinks_avg + compete_avg
sust <- sust1518 ~ sust95 + lgdpnrPC9095 + lrentPC9095 + dem_avg +
  parlinks_avg + compete_avg
sur_sys <- list(incl = incl, sust = sust)

sur_res1 <- systemfit(sur_sys, data = subset(exis, !is.na(sust1518)),
                      method = "SUR")
summary(sur_res1)

sur_res_lm <- systemfit(sur_sys, data = subset(exis, !is.na(sust1518)),
                      method = "OLS")
lrtest(sur_res1, sur_res_lm)  # likelihood-ratio test (SUR vs. OLS)

